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Abstract: We determine the spatial (impact parameter) dependence of nuclear parton 
distribution functions (nPDFs) using the ^-dependence of the spatially independent (av- 
eraged) global fits EPS09 and EKS98. We work under the assumption that the spatial 
dependence can be formulated as a power series of the nuclear thickness functions Ta- 
To reproduce the ^4-dependence over the entire x range we need terms up to [Ta] 4 - As 
an outcome, we release two sets, EPS09s (LO, NLO, error sets) and EKS98s, of spatially 
dependent nPDFs for public use. We also discuss the implementation of these into the 
existing calculations. With our results, the centrality dependence of nuclear hard-process 
observables can be studied consistently with the globally fitted nPDFs for the first time. 
As an application, we first calculate the LO nuclear modification factor R^a f° r primary 
partonic-jet production in different centrality classes in Au+Au collisions at RHIC and 
Pb+Pb collisions at LHC. Also the corresponding central-to-peripheral ratios R l (jp are 
studied. We also calculate the LO and NLO nuclear modification factors for single inclu- 
sive neutral pion production, -R^Au' at mid- and forward rapidities in different centrality 
classes in d+Au collisions at RHIC. In particular, we show that our results are compatible 
with the PHENIX mid-rapidity data within the overall normalization uncertainties given 
by the experiment. Finally, we show our predictions for the corresponding modifications 
i?pP b in the forthcoming p+Pb collisions at LHC. 
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1 Introduction 

In a high-energy hadronic or nuclear collision of particles A and B the inclusive cross 
sections for hard processes involving a large interaction scale Q 2 3> Aq CD can be computed 
using the QCD collinear factorization theorem [1, 2], 

da AB^k+x = £ f A (Q 2 ) ® ff{Q 2 ) ® da l ^ k+x ' + 0(1/Q 2 ), (1.1) 

where d<r represents the perturbatively computable partonic pieces (cross sections in lowest 
order), and iff) is the parton distribution function (PDFs) for a given parton flavor i in 
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the colliding particle A (and correspondingly for the flavor j in B). The PDFs are universal, 
process-independent functions of nonperturbative origin, whose evolution in the scale Q 2 
can, however, be obtained from the DGLAP equations [3-6] derived from perturbative 
QCD. 

A precise knowledge of the universal PDFs is thus vital for interpreting any hard- 
process results at the present colliders BNL-RHIC and CERN-LHC. This holds as well for 
proton-proton collisions as for proton- nucleus and nucleus- nucleus collisions. To determine 
the nonperturbative input in the PDFs, one has developed global analyses which exploit 
a multitude of experimental hard-process data and the DGLAP evolution. Excellent fits 
for the free proton PDFs have been obtained, and sets like CT10 [7], MSTW [8], and 
NNPDF2.0 [9] are nowadays available. 

It is well known that the PDFs of nucleons bound to a nucleus, the nuclear PDFs 
(nPDFs), are modified relative to the free-nucleon PDFs. Analogously to the free-proton 
case, global DGLAP analyses have been developed also for the nPDFs. The further com- 
plication with these is that in addition to the usual x and Q 2 dependences also the nuclear 
mass-number (^4) dependence of the PDFs needs to be dealt with. The global nPDF fits 
have so far resulted in leading-order (LO) nPDF sets EKS98 [10], HKM [11] and HKN04 
[12], and next-to-leading order (NLO) sets nDS [13], HKN07 [14], EPS09 [15], nCTEQ 
[16, 17], and DSZS [18]. Importantly, and similarly to the free proton case, with the error 
sets of EPS09 (and similar sets in DSZS), one can nowadays quantify how the uncertainties 
remaining in the nPDFs, illustrated in Fig. 1, are transmitted to the nuclear hard-process 
cross-sections. 

The global analyses mentioned above have all considered only the spatially averaged 
nPDFs, probed in minimum-bias nuclear collisions with no cuts on the collision centrality 
(impact parameter). In particular, as the modest amount of available nuclear hard-process 
data severely limits the number of possible fit parameters, it has so far been impossible to 
embed the spatial dependence, or the impact-parameter dependence, of the nPDFs directly 
into the global analysis. An obvious drawback with the globally analysed nPDFs then is 
that it has not been possible to consistently compute nuclear hard-process cross-sections 
in different centrality classes. 

The purpose of this study is to consider this problem by pinning down the spatial 
dependence of the nPDFs, i.e. flic dependence of the nuclear modifications of the PDFs 
on the nucleon's position in a nucleus. We do this in a manner which is for the first time 
fully consistent with the nPDFs from a global analysis. Earlier attempts to this direction, 
lacking however such a consistency, can be found in Refs. [19-22] . A further motivation for 
the current study is the Gribov-Glauber modeling of nuclear shadowing, reviewed lately in 
Ref. [23], whose output nPDFs are not a result of a global analysis like EPS09 but which 
have so far been the only ones where the spatial dependence arises in a self-consistent 
manner from modeling the physics origin of the nuclear effects. On the experimental side, 
the current study is inspired by e.g. the measurements of single hadron production [24-30] 
and J/^> production [31, 32] in different centrality classes in d+Au collisions at RHIC, as 
well as by the hard-process measurements in the forthcoming p+Pb collisions at the LHC. 
Also the theoretical modeling of the J/^f production discussed recently in Ref. [33] has 
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Figure 1. The nuclear modifications and their uncertainties in a lead nucleus (A = 208) for 
different parton flavors from EPS09NLO at the EPS09 initial scale Qq = f .69 GcV 2 (upper panel), 
and from EPS09LO and EKS98 at the EKS98 initial scale Q\ = 2.25 GeV 2 (lower panel). 



motivated our study. 

Our basic idea for uncovering the spatial dependence in the EKS98 and EPS09 nPDFs 
is straightforward: We first introduce the spatial dependence of the nuclear modification 
to the nPDF of each parton type i in each nucleus A at each x and Q 2 in terms of a power 
series of the standard nuclear thickness functions Ta- Then, we determine the coefficients 
of each power of Ta by exploiting the A-dependence of the EPS09 and EKS98 nPDFs 
(these sets, through the global fits, represent the experimental data here). As an output, 
we provide the numerical routines named EPS09s (LO and NLO as well as error sets for 
both) and EKS98s for computing the spatially dependent nPDFs which - simultaneously 
for all nuclei considered - normalize to the corresponding spatially independent EPS09 and 
EKS98 nPDFs. These new sets will be downloadable at the link [34]. 

As concrete examples of how to easily implement our spatially dependent nPDFs and 
the nuclear collision geometry in the computation of nuclear hard-process cross-sections 
in different centrality bins, we first discuss the centrality dependence of the LO nuclear 
modification ratios J^aa^Pt) °^ P r i mar y partonic-jet production in Au+Au collisions at 
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RHIC and Pb+Pb collisions at the LHC. We also study the nuclear modification factors of 
inclusive 7r° production, R^ Au , in d+Au collisions at RHIC and R^p h in p+Pb collisions at 
the LHC, both at mid- and forward rapidities, and considering both the NLO and LO cases. 
For -Rja u we also make, to our knowledge, a first comparison with the PHENIX centrality 
dependent data [26] where the overall normalization errors of the data are accounted for 
in detail. Due to the planned p+Pb program at the LHC, the ratio R^p^ipr) for single 
hadron production has been of growing interest recently [35-40], and we will show also 
here how interesting and useful this ratio would be from the point of view of constraining 
the nPDFs further. 

The paper is organized as follows: In Sec. 2 we define the model framework and explain 
the fitting procedure. In Sec. 3 we show the results for the spatially dependent nuclear 
modifications of PDFs. Also a comparison with selected other works is presented here. 
Applications of our results are discussed in Sec. 4. For clarity, a summary of the standard 
elements used in the applications here, the formulation of the nuclear collision geometry, 
different overlap functions and the optical Glauber model, is given in the Appendix A. 

2 The Analysis Framework 

2.1 Definitions of the Nuclear Modifications 

First we need to define how we introduce the spatial dependence to the nPDFs in terms 
of the hard-process cross-sections. Let us start with the usual spatially averaged nPDFs. 
The number distribution of an observable k produced in a collision of nuclei A and B at 
an impact parameter b is given by 

dN AB ^ k+x (h) = T AB (h)da AB ^ k+x , (2.1) 

where TabO*) is the standard nuclear overlap function normalized to AB (cf. Eq. (A. 20) 
in App. A. 2, see the nuclear collision geometry in Fig. 20), and da AB ^ k+x is the b- 
independent inclusive hard cross-section of Eq. (1.1) containing the nPDFs and pertur- 
bative pieces. The spatially averaged nPDFs in a nucleus A with Z protons and A — Z 
neutrons are now given by 

= jf? /A (x,Q 2 ) + ^f? /A (x,Q% (2.2) 

where the nPDFs of a bound neutron, f^ A , may be (approximately) obtained from those 
of the bound proton, ff , by using the isospin symmetry (see [15]). As in EKS98 and 
EPS09, we define the nPDF for each parton flavor in terms of the spatially averaged nuclear 
modification R A (x,Q 2 ) and the corresponding free proton PDF ff(x,Q 2 ), 

f? /A (x,Q 2 ) = R A (x,Q 2 )ff(x,Q 2 ). (2.3) 

To lighten the notations, we express the nPDFs in Eq. (2.2) as 

f t A (x,Q 2 ) = ^i?P(x,Q 2 )/f (x,Q 2 ), (2.4) 
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where the sum runs over all the nucleons N = 1, . . . , A. 

Decomposing the Tab into the standard nuclear thickness functions (cf. Eq. (A. 20)), 
and using Eq. (1.1) we may write 



dN AB ^ x (b) = £ -L £ fd 2 s 1 T A (s 1 )^ M (x 1 ,Q 2 )/^(x 1 ,Q 2 ) 1 



/" d 2 s 2 T B (s 2 ) Rf B/B (x 2 ,Q 2 ) ff B (x 2 ,Q 2 ) ® d<r^ fc+x '5(s 2 - si - b). 



(2.5) 



From this, we see that a suitable definition of the spatially dependent nuclear modification 
r A (x,Q 2 ,s) for the PDF of parton flavor i (per nucleon) is 

R A (x,Q 2 ) = 1 / d 2 sT A (s) ? f (x,Q 2 ,s), (2.6) 

where the thickness function is normalized to j4 and where the case of no nuclear effects 
corresponds to R A = rf = 1. Using these definitions, we can now generalize Eq. (2.5) to 
include the spatially dependent nuclear modifications, 



dN AB ^ x (b)=J2^B £ [d 2 s 1 T A (s 1 )rf(x 1 ,Q 2 ,s 1 )f-(x 1 ,Q 

i,j,X' N A ,N B J 

J d 2 S2 TB(s 2 )rf(x 2 ,Q 2 ,s 2 )ff B (x 2 ,Q 2 )®da^ k+x '5(s 2 - Sl - b). 



(2.7) 



As a consistency check, we note that the definition in Eq. (2.6) guarantees that the 
minimum-bias cross sections, which are obtained by integrating Eq. (2.7) over the whole 
b space, become simply AB times the hard cross-section computed with the spatially 
averaged nPDFs, 

^mT k+X = [ d 2 bdN AB ^ k+x (b) = AB f\ A (x,Q 2 )®ff(x,Q 2 )®da^ k+x '. (2.8) 

The key assumption in the present analysis is that the spatial dependence of r A (x, Q 2 , s) 
is a function of the nuclear thickness ?a(s). The motivation for this comes mainly from 
the shadowing region at small x, where the partons of sufficiently small values of x may 
interact with partons from any other nucleon near enough in the transverse direction. Also 
in the Gribov-Glauber modeling [23] of the initial state nPDFs the nuclear effects become 
essentially functions of Ta- The functional form we choose to use and test here is a simple 
power series of the thickness functions, 



r A 
' i 



(x, Q 2 , s) = 1 + ]T dj(x, Q 2 ) [T A (s)Y ■ (2-9) 

3=1 



Here we would like to emphasize the following points: First, all the A dependence 
is now in the thickness functions which are fully known, and all the coefficients dj(x,Q 2 ) 
which will be our fit parameters, depend on x and Q 2 but not on A. Second, the power 
series of the form 1 + . . . also fixes by construction that r A (x,Q 2 , s) —> 1 when |s| — > oo, 
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which means that the nucleons at the very edge of the nucleus are essentially regarded as 
free nucleons. Third, what is known from the EKS98 and EPS09 -types of analyses, are 
only the spatially averaged nuclear modifications and their A systematics, i.e. T^-weighted 
integrals of Eq. (2.9) over s for each nucleus. Fourth, since the EKS98 and EPS09 global 
analyses have not been constructed to reproduce any specific theoretically motivated A 
dependence of the nPDFs, we can test the validity of the assumption of Eq. (2.9), as well 
as the number of terms needed, only a posteriori. 

Using the definitions above, we can see why the simplest 1-parameter approach with 
n = 1 in Eq. (2.9) (which is used e.g. in [19-22] as well as in e.g. the HIJING event 
generator [41]) is not fully consistent with the observed A systematics of the nuclear data. 
In this case, rf(x,Q 2 ,s) = 1 + c*(x, Q 2 )Tk(s), and from the definition in Eq. (2.6), one 
obtains c^Q 2 ) = [Rf(x,Q 2 ) - 1]A/T AA (0), where Rf is given by the globally analysed 
nPDFs (i.e. nuclear data). The problem then is that the coefficient c l (x,Q 2 ) may depend 
in fact quite strongly on A, which indicates that the simplest assumption of terminating the 
power series at the first nontrivial term does not correctly capture the spatial dependence 
of the measured nuclear structure functions. This redundant A dependence is illustrated 
in Fig. 2 for gluons in a lead nucleus at x = 0.01 at the initial scales of the sets EKS98, 
EPS09LO1 and EPS09NLO1. We can see that especially for the NLO set the problem is 
more serious. One of the driving motivations for the present study is to solve the problem 
of recovering the A systematics in the spatially dependent nPDFs. 
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Figure 2. The problematic A dependence of the parameter c 9 (x,Q 2 ) ~ [R^(x,Q 2 ) — 1]A/Taa(0) 
for EPS09NLO1 and EPS09LO1 (EKS98) gluons at x = 0.01 and Q 2 = 1.69 (2.25) GeV 2 in the 1- 
paramctcr approach where one includes only the first nontrivial term in the power series in Eq. (2.9). 

2.2 Fitting Procedure 

To extract the ^4-independent coefficients c*(x,Q 2 ), we need to introduce a fitting proce- 
dure, where we utilize the definition (2.6) and the A dependence of the EKS98 and EPS09 
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nuclear modifications at different values of x and Q 2 for each parton flavor i. To reproduce 
the A systematics in the spatially independent nuclear modifications with the power-series 
ansatz of Eq. (2.9), we minimize the x 2 defined as 



where the spatially averaged modifications Rf(x,Q 2 ) from EKS98 and EPS09 now rep- 
resent the "experimental" data. The weight factors W^(x,Q 2 ) are artificial errors which 
control the quality of the fit and which are set by hand. Our numerical observation is that 
for good fits we need 4th-order polynomials in Ta, i.e. n = 4 in Eq. (2.9). Furthermore, best 
fits were obtained with the weig ht Wf(x, Q 2 ) = Rf{x, Q 2 ) - 1 for EKS98 (this corresponds 
to fitting the deviations from unity within a constant relative error) and Wf{x,Q 2 ) = 1 
for EPS09 (corresponds to fitting the modifications within a constant error). 

By construction, both EKS98 and EPS09 give no nuclear modifications for deuterium. 
This cannot be reproduced with the functional form we selected for rf(x,Q 2 ), and we do 
not expect the fit form of Eq. (2.9) work for the smallest values of A, either. Consequently, 
we exclude the nuclei A < 16 from the fit. Thus for EKS98 the sum runs over A = 
16, 20, ... , 300 (i.e. emphasising the large nuclei) and for EPS09 we use all the A > 16 
values for which these sets are currently available. 

3 Results 

3.1 Quality of the fit 

First, we demonstrate that our fit framework manages to reproduce the spatially averaged 
nuclear modifications and especially their A dependence indeed very well. Figure 3 shows 
the obtained spatially dependent gluon modifications integrated over the transverse plane 
according to Eq. (2.6), and the corresponding input modifications at different fixed values of 
Q 2 from the NLO set EPS09NLO1 (left panel), and from the LO sets EPS09LO1 and EKS98 
(right panel), for a lead nucleus. In what follows, we refer to these cases as "EPS09sNLOl" , 
"EPS09sLOl" and "EKS98s" where "s" is for "spatial" and "1" for the central sets. As 
seen in the figure, the match with the input and output distributions is very good; for all 
parton flavors and the nuclei included in our fits it is within 2 % at x < 0.75 for EPS09NLO, 
1 % at x < 0.85 for EPS09LO, and 0.2 % at x < 0.95 for EKS98. Importantly, the key- 
feature here, the A dependence of EPS09 and EKS98, is similarly well reproduced, as is 
demonstrated by Fig. 4 below. 

Recall also that in the EPS09 global analysis in addition to the best fit there are also 
30 error sets, which enables one to compute how the uncertainties of the nPDFs propagate 
into physical observables. The above fitting and determination of the spatial dependence 
are done also for each of these error sets, both in LO and in NLO, and the fit quality is 
similar as in Figs. 3 and 4. Thus the error propagation calculations (as instructed in [15]) 
for centrality-dependent nuclear hard cross-sections can now be done as before, using the 
EPS09s sets. 




(2.10) 
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Figure 3. Left: The spatially averaged nuclear modification Rf(x, Q 2 ) for a lead nucleus (A = 208) 
from the NLO set EPS09NLO1 (dotted lines) and from the EPS09sNLOl spatial fit presented here 
(solid lines) at four different scales. Right: The same with the LO sets EKS98 and EPS09LO1 
(dotted) and with the spatial fits EKS98s (dashed) for three different scales and EPS09sLOl (solid) 
for four different scales. 




Figure 4. Left: The A dependence of the spatially averaged nuclear modification R^(x,Q 2 ) at 
fixed values x = 0.001 and Q 2 = 1.69 GeV 2 from the sets EPS09NLO1 (crosses) and EPS09LO1 
(pluses) and from the corresponding spatial fits EPS09sNLOl (solid green line) and EPS09sLOl 
(solid blue line). Right: The same but with the LO set EKS98 (circles) and the corresponding 
spatial fit EKS98s (solid red) at Q 2 = 2.25 GeV 2 . The small nuclei shown with gray markers in 
both panels were not used in our spatial fits. 



3.2 Spatial Dependence 

After the consistency checks above, let us next discuss the spatial dependence obtained 
for the nuclear modifications of the PDFs. In Fig. 5 we present the nuclear modification 
r^ b (x, Q 2 , s) at the initial scale Q 2 = 1.69 GeV 2 as a function of x and s, as obtained from 
the fitting to the sets EPS09 NLO and LO, as well as the LO set EKS98. The three main 
observations are 

• The overall x-shape of the nuclear modification away from the edge of the nucleus, 
at |s| < Ra, is similar as in the input distribution. This confirms that our fit does 
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not generate any unwanted extra curvature. 



• In the center of the nucleus, |s| ~ 0, the nuclear modification is only slightly larger 
than the input average modification. This also confirms the earlier similar findings 
in [19, 20, 22]. 

• The nuclear modification dies out as expected, by construction, when |s| > Ra- This 
feature arises from the vanishing 7a(s) at the edge of the nucleus. 

The observations for the spatial dependence of the sea and valence quarks nuclear modifi- 
cations are the same. Examples of these can be found in App. B. 

3.3 Comparison with other approaches 

Next, we compare our EPS09s and EKS98s fits with the 1-parameter approach described in 
the end of Sec. 2.1. [20, 22]. 1 This model has been used to study the centrality dependence 
of the J suppression e.g. in Refs. [31, 33, 42, 43] 2 and inclusive hadron production in 
d+Au collisions at RHIC in Ref. [22]. We also compare with the leading- twist formulation 
[23, 44] of nuclear shadowing which is based on the generalization of Gribov-Glauber theory, 
QCD factorization and diffractive PDFs measured at HERA. For the spatially averaged 
nuclear modifications, this model typically predicts a stronger smallest-x shadowing than 
what is implemented in the parametrizations of EKS98 and EPS09 (see e.g. Ref. [23]). For 
the comparison, we consider the FGS10_L set [23, 45], and choose the value of x not too 
small, so that the spatially averaged FGS10_L nuclear gluon modification is close to that 
in EPS09 or EKS98. 

In Fig. 6 we plot the nuclear modification for gluons at fixed values of x and scale Q 2 = 
4 GeV 2 for A = 208 as a function of |s| from our EPS09sNLOl, EPS09sLOl and EKS98s 
fits, from the 1-parameter approach using the averaged sets EPS09NLO1, EPS09LO1 and 
EKS98, and from FGS10_L. Although numerically the differences are not very large, we 
notice that while both the EPS09sNLO and EKS98s results are close to FGS10_L, the 
1-parameter approach leads to a too steep transverse profile for the modifications in all 
cases. 

4 Applications 

Next, we consider some concrete examples of computing the nuclear hard-process cross- 
sections in different centrality classes using the spatially dependent nPDFs. First, we 
discuss the centrality dependence of primary partonic-jet production in A+A collisions at 
RHIC and LHC. Then, we consider neutral pion production in d+Au collisions at RHIC 
and in p+Pb collisions at the LHC. 

1 In [20] the spatial dependence enters through the first nontrivial power of the nuclear density Pa{v) or 
the thickness function Ta(s). The latter scenario corresponds to what we refer to as " 1-parameter approach" 
here. 

2 In [33] one studies also other types of spatial dependences. 
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Figure 5. The spatially dependent modification of gluon distribution in a lead nucleus, 
r^ h (x,Q 2 ,s), from EPS09sNLOl (upper left), EPS09sLOl (upper right) and EKS98s(lower plot) 
as a function of x and s at the initial scale Q 2 = 1.69(2.25) GeV 2 of EPS09 (EKS98). For examples 
of the corresponding plots of other parton flavors, see App. B. 



4.1 Implementation of EKS98s and EPS09s 

For denning the centrality classes we use the optical Glauber model specified in App. A. 2. 
In this case, each centrality class corresponds simply to a certain impact-parameter interval 
|b| € [61,62]- The generic average number distribution of a hard-process observable k in 
this centrality class of an A+B collision is 

. v (, b2 d 2 bd^ B (b) 

diVlA , = ™: > (4.1) 



' A R 7 — 1 , 
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Figure 6. Comparison of the spatial dependence of the gluon modification in a lead nucleus, 
r^ h (x, Q 2 , s), between FGS10_L (short-dashed blue curves), 1-parameter approach (long-dashed 
green) and our spatial fits (solid red) EPS09sNLOl (upper left), EPS09sLOl (upper right) and 
EKS98s (lower plot). The scale Q 2 = 4 GeV 2 for all plots but the values of x have been chosen so 
that the spatially averaged R^ h (x,Q 2 ) (dotted horizontal red lines) approximately coincides with 
FGS10_L (dotted blue). 



where (b) = 1 — exp[— Tab(^)<^nn\ f rom Eq. (A. 23), and dN^ B (h) is obtained from 
Eq. (2.7). Using the expansion of rf(x, Q 2 ,s) in powers of Ta from Eq. (2.9), the integrals 
over the impact parameter for the spatially dependent parts can be conveniently separated 
from the spatially independent fit coefficients, free nucleon PDFs and pQCD parts as 
follows: 

4 i 

d 2 bdivi B (b)= y, n%(biM E jb £ <i( x uQ 2 )f? A ( x hQ 2 )® 

n,m=0 iJ,X' N A ,N B 

4n(x2,Q 2 )f? B (x 2 , Q 2 ) da^ k+x ' (4.2) 

where f^j A,Ns are the free nucleon PDFs, and we have defined c^(x,Q 2 ) = 1 and 

62 

TTB{hM) = I d 2 b J d 2 s [T A (s - b/2)] n+1 [T B (s + b/2)] m+1 . (4.3) 
61 
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From Eq. (4.2) we see the most straightforward implementation of the spatially de- 
pendent nPDFs. The purely geometric integrals, the ^^(61,62) in Eq. (4.3) for each pair 
of the powers n and m, can be computed independently of the kinematic variables x\, 
X2, and Q 2 and also independently of the parton flavors Thus, in total we have 25 
different geometric integrals to do (or 15 if A = B) but we need to do them only once. 
In comparison with the spatially averaged case, the fit parameters c l n (x, Q 2 ) thus play the 
role of the nuclear modifications R\(x,Q 2 ) for each of the 25 pairs n,m. To arrive at 
the final b-integrated result for the number distribution of k, we thus need to repeat the 
computation of the kinematic parts 25 times, each with different sets of the coefficient 
pairs {c4},{c4} and a different geometric weight T A l1 g(b 1 ,b 2 ). The EKS98s and EPS09s 
routines which we provide in [34], give in addition to the fit coefficients {d n (x, Q 2 )} also 
the thickness functions Ta(s) (used in the fits here) for the computation of Tj^g b 2 ), as 
well as the combination TA(s)rf(x, Q 2 , s) for other possible implementations. Note also 
that for the b integral in Eq. (4.3) the angular part is trivial, giving just 2ir. 

4.2 The Nuclear Modification Factors R 1 ]^ and R^f 

Let us now consider the centrality dependence of primary inclusive high-p^ parton pro- 
duction in A+A collisions at RHIC and LHC. Following the generic discussion above, we 
define the nuclear modification ratio Rj*£(pr) relative to the p+p case for each centrality 
class as 

/ d 2N AA \ t H 2 /V ljet (Vl 

i«( W ,K».,w=— ^ — rf^-— — *^*r. («) 

where (NS^r )b 1 & 2 is the average number of binary collisions in this centrality class given by 
Eq. (A.27) and 

^inei ^ khs inelastic nucleon-nucleon cross section. Apart from the (small) 
isospin effect, this ratio yields unity if there are no nuclear effects in the nPDFs. Thus, for 
peripheral enough centrality bins, this ratio should approach unity. For the details of the 
partonic cross sections, bookkeeping and kinematics, we refer to [46]. 

The nuclear mofication factor in the minimum-bias collisions is obtained from above 
by setting b\ = and 62 - > 00, in which case we have 

1 2 ljet r, ]j c t 

'oljet/ \\ 1 a a AA,MB /dV p J p 



"AA^^yJ/ A 2 dprdy I dprdy 

where ^Jambi wri ich contains only the spatially averaged nPDFs, is obtained from 
Eq. (2.8) by setting B = A, and the p+p baseline dopp* from the same equation by 
setting A = B = p. 
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In addition to the centrality dependence of Rl^i we are interested in the central-to- 
peripheral ratios, defined as 




where C and P refer to the central and peripheral bins, correspondingly. The advantage 
of this ratio (in the experiments) is that the information of the proton-proton baseline 
is not required. In particular, we would like to see exactly how much R C p differs from 
the modification RFa which is computed with the spatially averaged nPDFs. We perform 
these example-calculations for both RHIC and LHC but for simplcity only to LO pQCD, 
since without jet quenching these ratios do not directly correspond to observables. They 
illustrate, however, the points we wish to make with the spatially dependent nPDFs, and 
also serve as (LO) pQCD baselines for the observed suppression of high-p^ particles. 

The two different centrality classes we consider here for Au+Au collisions at RHIC 
and Pb+Pb collisions at the LHC, are the central 0-5% and peripheral 60-80% bins. The 
Glauber model input and the resulting impact parameter intervals and average numbers 
of binary collisions in these centrality classes are summarized in Table 1 . 



Table 1. The centrality classes as impact parameter intervals, and average number of binary 
collisions from the optical Glauber model in A+A collisions for RHIC and LHC. 



/ NN 

[GeV] [mb] 


Central = 0-5 % 
h [fm] b 2 [fm] (N bin ) 


Peripheral = 60-80 % 
h [fm] b 2 [fm] (N bin ) 


Au+Au 200 42 
Pb+Pb 2760 64 


0.0 3.355 1083 
0.0 3.478 1771 


11.62 13.42 15.10 
12.05 13.91 19.08 



In Fig. 7 we plot the ratio Rl% {pt, y = 0) for central, peripheral and minimum-bias 
collisions, as well as R^plpx) in Au+Au collisions at RHIC. Figure 8 shows the same 
quantities for the LHC Pb+Pb case. The central and peripheral RJa and R^p have been 
obtained with the spatially dependent nPDFs EPS09sLOl (left) and EKS98s (right), and 
the average (R>]Ia) m minimum bias collisions with the spatially independent EPS09 and 
EKS98 nuclear modifications. For the free proton PDFs we have used CTEQ6.1L [47]. 
The renormalization scale [i and factorization scale Q has been set to be the transverse 
momentum, pr, of the parton. 

The main observations from the figures are: (i) The central Rj^ is quite close to the 
average RJai which is expected since the nuclear modifications at small s are close to the 
average modifications, (ii) The peripheral RFa ls clearly not unity but there appear almost 
10% antishadowing effects at mid-pT at RHIC and even more than 20% shadowing effects 
at small pt at the LHC, and up to 10% EMC effects at large px both at RHIC and LHC. 
(Hi) Consequently, the ratio Rqp differs significantly from the average Rj^- The results 
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Figure 7. The LO nuclear modification Rjj± as a function of partonic transverse momentum for 
central (red long-dashed), peripheral (blue dashed) and minimum-bias (green dot dashed) collisions, 
and Rqp (solid magenta) for Au+Au collisions at ^s^n = 200 GcV and y = using EPS09sLOl 
(left panel) and EKS98s (right panel) . 




2 5 10 20 50 100 200 500 2 5 10 20 50 100 200 500 

p T [GeV/c] p T [GeV/c] 

Figure 8. The same as Fig. 7 but for Pb+Pb collisions at ^snn = 2.76 TeV 



suggest that in a precision theory-analysis of the centrality dependence of jet quenching, 
one needs to account also for the spatial dependence of nPDFs. Finally, regarding the 
differences between the different nPDF sets applied here, we observe in Figs. 7 and 8 how 
the stronger shadowing in the EPS09 case (cf. Fig. 3) translates into steeper px slopes of 
RFa at small px than in the EKS98 case. 

4.3 Centrality dependence of R^ u {pt) at RHIC — comparison with data 

While the above ratios Rjj^ and R]jp mainly serve as theoretical pQCD baselines for jet 
quenching studies, it is important to test our spatially-dependent nPDF framework against 
some measured centrality-dependent observables. To avoid the complications of hot QCD 
matter modeling, we turn to the highest-energy d+Au collisions at RHIC and p+Pb at the 
LHC. For our purposes a promising published data set is the nuclear modification factor 
RdAuiPT) f° r single inclusive neutral-pion production at mid-rapidity |?7| < 0.35, measured 
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by PHENIX [26] at different centrality classes in d+Au collisions at y / s/v r /v = 200 GeV. 
Since the minimum-bias data precisely from this data set was among the constraints in the 
global EPS09 fits, it is now very interesting to study, for the first time consistently with 
EPS09, how well we can reproduce the measured centrality dependence of this ratio. 

Analogously with Eq. (4.4), we define the centrality-dependent nuclear modification 
factors as 

R&^y-MM^ r§^ = — ^f^' M 

where the number distribution now involves a further folding over the fragmentation func- 
tions, 

dJV3>) = dN LQ>) ® D n0/k (z, Ql) (4.8) 
k 

and where &N^ A (h) is obtained from Eq. (2.7) by setting A = d and B = A, and (as we do 
not assign any nuclear effects for the deuterium PDFs) also rf = 1. For obtaining a realistic 
thickness function for deuterium, we use the Hulthen wavefunction formulation [48] given 
in App. A. 1.2. The impact parameter ranges and average numbers of binary collisions at 
the corresponding centrality classes are obtained again from the optical Glauber model. 

4.3.1 Minimum-bias R^^ipr) 

Setting the spatial integrals in Eq. (4.7) over the whole impact-parameter space gives again 
the minimum- bias ratios, 

r ]2 -it j 2 7r° 

R ^PT,y)) = ^^^-/^, (4.9) 

where da^° AMB = Ylk da dAMB ® /k{Z)Q%) again contains only the spatially averaged 
nPDFs in da k AA 

mb which is obtained from Eq. (2.8) by setting A — d, B — A. As noted 
earlier, in the EKS98 and EPS09 frameworks there are no nuclear modifications to the 
deuterium PDFs. The p+p baseline diTp" is computed correspondingly, but without any 
nuclear effects. 

Figure 9 shows the PHENIX data [26] and our NLO (left) and LO (right) results 
for the nuclear modification factor (-RJa. u (pT) D = 0)) in minimum-bias collisions. For 
the NLO calculation with EPS09sNLOl (equivalently one may use EPS09NLO1, since 
the spatial dependence is here irrelevant) we used the NLO fragmentation functions from 
KKP [49] 3 , AKK [50] and fDSS [51]. For the free proton PDFs, we use CTEQ6M [47]. 
Correspondingly, the LO case was computed with EPS09sLOl and EKS98s, using the KKP 
and fDSS LO fragmentation functions and CTEQ6.1L PDFs [47]. The renormalization 
scale fi, factorization scale Q and fragmentation scale \xf are all fixed to pr, the transverse 
momentum of the produced hadron. For details of the LO calculation, we again refer to 
Ref. [46], while the NLO computation was performed by using the INCNLO code [52, 53]. 



The KKP set was also used in the EPS09 global analysis. 
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Figure 9. The nuclear modification factor R^^ipr) for ^/s/vat = 200 GcV at y = for minimum 
bias collisions. Calculations are for the NLO pQCD using EPS09NLO1 with three different frag- 
mentation functions (left), and LO using EPS09LO1 and EKS98 with two different fragmentation 
functions (right). The blue bands are computed using the EPS09 error sets and fDSS fragmentation 
functions. The experimental PHENIX data [26] are shown by markers, and their error bars (boxes) 
stand for the point-to-point statistical (systematic) errors. Notice that the data points and their 
errors have been multiplied by a factor 1.039 (1.050) for the NLO (LO) case, which is well within 
the 9.7 % overall normalization error quoted by the experiment (see text for details). 



From Fig. 9, we notice the following: (i) The EPS09 uncertainty bands for the NLO 
results are slightly smaller than in the LO case, reflecting the fact that the EPS09NLO glu- 
ons are somewhat better constrained in the antishadowing region than those of EPS09LO. 
(ii) In the small px region there is a difference in the pt slopes between the EKS98 and 
EPS09LO1 results. This is caused by the weaker shadowing in EKS98. However, also the 
EKS98 results remain within the EPS09 error bars. (Hi) The uncertainty caused by the 
differences in the fragmentation functions remains conveniently small in all cases. 

Regarding the data comparison in Fig. 9, we emphasize the following important point: 
In addition to the the statistical uncertainties (error bars) and point-to-point systematic 
errors (boxes), PHENIX quotes a 9.7 % overall uncertainty which originates from the p+p 
reference and which is not included in the statistical error bars shown. Consequently, 
allowing for a shift of the data points and their errors by less than 9.7 % and requiring the 
best possible overall fit to the data (using the fDSS FFs and by minimizing the x 2 with the 
point-to-point statistical and systematic errors added in quadrature), we have multiplied 
the data by a factor 1.039 (NLO) and 1.050 (LO). Such a few-percent shift is well within 
the uncertainty given by the experiment. As already noticed in the EPS09 analysis [15], 
the resulting agreement with the data is quite good, both in LO and in NLO. 

Figure 10 shows the ratio {R^^iPTiU = 3)) at the forward region in minimum- bias 
collisions. Again, we show both the NLO (left) and LO (right) results with the same set-up 
as in Fig. 9 above. Now the differences between the fragmentation functions start to be 
visible, as one is probing their larger- z tails where the uncertainties are larger: for the same 
pion pt, the differences in the large- z fragmentation functions map to different values of x 
in the nPDFs. We also notice that the LO calculation gives a stronger small-pr suppression 
than the NLO case, which is partly due to the different pace of the scale evolution with the 
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NLO and LO nPDFs (cf. Fig. 3) and partly because the NLO computation probes slightly 
higher values of x than the LO case. Like in Fig. 9, the EPS09 error band is smaller for 
NLO than for LO. 4 



min. bias 
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Figure 10. The same as Fig. 9 but for neutral pions at a forward rapidity y = 3. 



4.3.2 Centrality dependent R^ Au {pt) 

Let us then look at the centrality dependence of the ratio R^° AvL - Our NLO and LO results 
for Rf An in the centrality bins - 20%, 20 - 40%, 40 - 60% and 60 - 88% are plotted 
in Figs. 11 and 12, correspondingly, together with the PHENIX data. Table 2 lists the 
impact parameter ranges and average number of binary collisions for each centrality class, 
obtained from the optical Glauber model with cr^[ = 42 mb. 

Table 2. The centrality classes as impact parameter intervals, and average number of binary 
collisions from optical Glauber model for d+Au collisions at tJsnn = 200 GeV using = 42 mb. 







bi [fm] 


b 2 [fm] 


{Nun) 





-20% 


0.0 


3.798 


15.57 


20 


-40% 


3.798 


5.371 


10.95 


40 


-60% 


5.371 


6.583 


6.013 


60 


-88% 


6.583 


8.336 


2.353 



Again, it is important to consider the different overall normalization errors in the 
experimental data. For the centrality-dependent ratios plotted in Figs. 11 and 12 there still 
is the 9.7 % overall systematic error due to the p+p baseline discussed above. In addition 
to this, an overall normalization error of 6.6-9.6 % arising from the determination of the 
average number of binary collisions, is quoted separately for each centrality bin. Following 
again the same procedure as for Fig. 9, we multiply the data and their point-to-point errors 

4 We note that there exist experimental data for R^au a t larger rapidities [54, 55] which suggest a more 
substantial small-pr suppression than what the EPS09 and EKS98 predictions could accommodate. This 
deviation calls for a more detailed investigation which, however, is clearly beyond the scope of the present 
paper. 
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by a factor which minimizes the difference to our calculation. Even the largest upwards 
shift, 11.3 % for the centralmost bin in the LO case, is well within the acceptable total 
overall normalization error quoted by the experiment. Note also the systematic decrease 
of the multiplication factor from central to peripheral collisions, which we believe is due to 
the difference in the experimental and Glauber-model definitions of the centrality classes. 



2 4 6 8 10 12 14 16 2 4 6 8 10 12 14 16 
1.5 I 1 < 1 < 1 < 1 ■ 1 < 1 < 1 < 1 1 < 1 < 1 ■ 1 ■ 1 < 1 < 1 < 1 1.5 




p T [GeV/c] p T [GeV/c] 

Figure 11. The nuclear modification factor -RJauCpt) for -y/ijvw" = 200 GeV at y = for different 
centrality classes. Calculations are in NLO pQCD using EPS09sNLOl and three different fragmen- 
tation functions. The blue error bands are computed with the error sets EPS09sNLOx (x=2....,31) 
and fDSS, and the data are from PHENIX [26]. The set-up and labeling are the same as in the 
left panel of Fig. 9. Notice that the experimental data have been multiplied by a different factor 
in each panel, which all are well within the total overall normalization uncertainties given by the 
experiment (6.6, 6.7, 8.5, 9.6 % for the four centrality bins from the Glaubcrization and 9.7 % from 
the p+p baseline.) 

From Figs. 11 and 12 we observe that within the experimental and theoretical uncer- 
tainties our calculations are consistent with the measurements. Especially the centrality 
systematics obtained from our spatially-dependent nPDFs agrees quite well with the data: 
the nuclear modifications are strongest in the most central collisions and systematically 
weaken when going to more peripheral collisions. This is especially nicely reflected in the 
region 1.3 < px < 4 GeV, where the pt slopes (which are not affected by the overall mul- 
tiplications) become steeper towards more central collisions. We also see that, like in the 
minimum-bias case, the EPS09 error bands are slighty smaller for the NLO than for the 
LO case, and that the uncertainties arising from the fragmentation functions remain small. 

In Figs. 13 and 14 we plot also our NLO and LO results for -RJ^u a * V SNN = ^00 GeV 
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p T [GeV/c] p T [GeV/c] 

Figure 12. The same as in Fig. 11 but with LO calculations using EPS09sLOl and EKS98s and for 
two different fragmentation functions, and the blue error band is computed with the EPS09sLOx 
(x=2,...,31) and fDSS. Notice again the overall multiplicative factors for the experimental data. 

at a forward rapidity, y = 3, in the different centrality bins. In the forward region the nu- 
clear modifications are larger since we now are probing smaller x values in the nPDFs than 
in the mid-rapidity region. Like in the minimum-bias case, we notice that the difference 
between the fragmentation function sets we use, becomes noticeable in the forward region. 
Again the small-pT suppression is stronger and nPDF-orginating uncertainties are larger 
for the LO case. 

4.4 Predictions for p+Pb collisions at LHC 

In the heavy-ion program of the LHC at CERN, there are now plans to collide protons 
with lead nuclei. Such collisions would be very useful for testing the QCD factorization 
and the universality of nPDFs, as well as for constraining the nuclear PDF modifications 
further especially at small values of x. Also the centrality dependence of nPDFs could be 
examined in these collisions via inclusive hadron production, similarly to the RHIC d+Au 
collisions discussed above but without the theoretical uncertainties arising from modeling 
the deuterium geometry. Thus, it is interesting to see what are the predictions from our 
spatially dependent nPDFs for these collisions. 

In Fig. 15 we plot our EPS09sNLO results for the nuclear modification factor Rpp^ipr) 
for neutral pion production in p+Pb collisions at y/s^N = 5.0 TeV at y = in four different 
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Figure 13. The nuclear modification factor R^ n (pr) for ^sjvn = 200 GeV at y — 3 in different 
centrality classes. The computation is done in NLO pQCD using EPS09sNLOf and three different 
fragmentation functions. The error bands are computed with the EPS09sNLOx (x=2,...,31) and 
fDSS. 



centrality classes. 5 We use again the KKP, AKK and fDSS fragmentation functions here. 
The uncertainty bands arising from EPS09sNLO are computed using fDSS. The inelastic 
cross section er^[ = 70 mb for this y / J/vjv is obtained from Fig. 5 of Ref . [56] . This leads to 
the impact parameter values and the average number of binary collisions for each centrality 
class given in Table 3. For the projectile proton, we have not assumed any spatial size, so 
that relative to the deuterium case above, in the collision geometry we replace the thickness 
function Td(s) by S(s) and the overlap function TdA(b) by the thickness function Tpb(b). 

As can be seen from Fig. 15, the nuclear modifications are strongest in the small-py 
region in all centrality classes. To see the behaviour of -Rp Pb in this region more clearly, we 
plot the results also in logarithmic scale in Fig. 16. We again observe the general behavior 
which follows from the spatial dependence of the nPDFs: the nuclear modifications are 
stronger in the central collisions and weaker in the peripheral collisions. We also notice 
that the three fragmentation function sets yield almost identical results. 

Figure 17 shows the corresponding ratio in minimum bias p+Pb collisions, computed 
both in NLO (left) and in LO (right). Like in the forward-rapidity case at RHIC, and for 

5 Very recently, the LHC moved up to collisions energies ^/s pp = 8 TeV, hence we take ^snn = 
V s pp V Z/A S3 5.0 TeV. Note also that y is the rapidity in the NN cms frame, i.e. we do not include 
the rapidity shift due to the antisymmetric collision. 
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Figure 14. The same as Fig. 13 but for LO pQCD using EPS09sLOl and EKS98s with two 
different fragmentation functions. The error bands are computed with the EPS09sLOx (x=2,...,31) 
and fDSS. 



Table 3. The centrality classes as impact parameter intervals, and average number of binary 
collisions from optical Glauber model for p+Pb collisions at y/s = 5.0 TeV, with cr^3 = 70 mb. 







h [fm] 


b 2 [fm] 


{Nun) 





-20% 


0.0 


3.471 


14.24 


20 


-40% 


3.471 


4.908 


11.41 


40 


-60% 


4.908 


6.012 


7.663 


60 


-80% 


6.012 


6.986 


3.680 



the same reasons, the EPS09NLO leads to a weaker small-p^ suppression than EPS09LO 
and EKS98, and the uncertainty band is clearly smaller for the NLO case. 

As a probe of nuclear gluons even deeper in the small- x shadowing region, we plot in 
Fig. 18 our LO results for -Rp Pb at a forward rapidity, y = 3, for the four centrality classes. 
Again the KKP and fDSS fragmentation functions are used, and we see that they yield very 
similar results. We should also point out that the EPS09sLO error band for the peripheral 
bin in Fig. 18 can be regarded as an underestimate in that it has been computed without 

6 For y — 3, we could not obtain reliable results with INCNLO at pr < 5 GeV for this ^/snn, hence only 
the LO results are shown here. 
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Figure 15. The nuclear modification factor -RpPb(PT) for \fs = 5.0 TeV at y = for four different 
centrality classes, computed in NLO pQCD using EPS09sNLOl and three different fragmentation 
functions. The error bands have been obtained with EPS09sNLOx (x=2,...,31) and fDSS. 



1 2 5 10 20 50 100 200 2 5 10 20 50 100 200 




1 2 5 10 20 50 100 200 2 5 10 20 50 100 200 

p T [GeV/c] p T [GeV/c] 



Figure 16. The same as Fig. 15 but in a logarithmic scale to emphasize the small-pr region. 
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Figure 17. Left: The nuclear modification factor R^p\ } {pt) for ^Jsjsn = 5.0 TeV at y = 0, com- 
puted in NLO pQCD using EPS09sNLOl (left panel) and three different fragmentation functions. 
The error band is computed using EPS09sNLOx (x=2,...,31) and fDSS. Right: The same but in 
LO pQCD with EKS98s and EPS09sLOl with two different fragmentation functions, and the error 
band is for EPS09sLOx (x=2,...,31) with fDSS. 



the error set EPS09sLO7. The reason for this is that the error set EPS09LO7 gives in fact 
antishadowing at smallest x for the lightest nuclei, and in the EPS09sLO this maps into an 
antishadowing near the edges of a large nucleus. This unphysical feature can be cured only 
by redoing the EPS09LO global fit with an improved ^.-dependence of the fit functions. In 
the meantime, we suggest that a physically more meaningful upper limit for the LO error 
band in the small- x region for the peripheral bin can thus be obtained without this LO 
error set. 

Finally, in Fig. 19 we show the minimum-bias -Rpp b at y = 3 for the NLO case at 
Pt > 5 GeV and for the LO case starting from px = 1-3 GeV. Note the linear (logarithmic) 
Pt scale on the left (right). Again we notice the weaker suppression and smaller error 
bands in the NLO case. Comparing the right panels of Figs. 19 and Fig. 17, we see that 
the smallest-pT suppressions are of similar magnitude. This is because the ratio -Rpp b in 
the small-pr region at the LHC probes already at y = the flat part of the shadowing 
assumed as an input in EPS09 (cf. Fig. 4). Hence, a measurement of -Rpp b both in the mid- 
and forward-rapidities can be expected to serve as a relevant constraint for the smallest-x 
shadowing region. 

5 Summary and Conclusions 

We have developed a framework to determine the spatial dependence of the nuclear modi- 
fications of PDFs in such a way that the outcome is consistent with the globally analysed 
EKS98 and EPS09 nPDFs which in turn are DGLAP-based fits to nuclear hard-process 
data. Both the LO and NLO cases have been considered, and with EPS09 the spatial 
dependence has been extracted also for all the 30 error sets. Correspondingly, we call the 
obtained spatially dependent nPDF sets EPS09s and EKS98s. 

The spatial dependence is introduced in terms of powers of the nuclear thickness func- 
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Figure 18. The nuclear modification factor -RpPbCPr) for ^/sjvn = 5.0 TeV at y = 3 for four dif- 
ferent centrality classes, computed in LO pQCD using EPS09sLOl and two different fragmentation 
functions. The error band is computed using EPS09sLOx (x=2,...,31) and fDSS. 




Figure 19. Same as Fig. 17 but for a forward rapidity y = 3; the NLO (LO) results are on the left 
(right). 



tions Ta(s). Regarding the power series rf(x,Q 2 ,s) = 1 + c )( x iQ 2 ) [^k( s )] 5 > we have 

shown that the 1-parameter approach (n = 1, used e.g. in [19, 22]) is not sufficient for 
reproducing A systematics in the nPDFs, and that we obtain a good overall agreement with 
the globally analysed averaged nPDFs when we include terms up to [T4] 4 . The outcome 
of the performed fits, the sets of coefficients {Cj(x,Q 2 )} for each parton flavor i at each x 
and Q 2 , are tabulated separately for each of the nPDF sets we considered. These tables 
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along with a routine for interpolation and computing the needed thickness functions are 
downloadable at [34]. 

As a concrete application of our framework, we calculated the nuclear modification 
factor R]Ia for LO primary partonic jet production at different centralities in Au+Au col- 
lisions at RHIC and in Pb+Pb at the LHC. We observed that while the central R]^l is quite 
close to the minimum-bias ratio (R>]Ia) an d the peripheral R]Ia differs fairly significantly 
from unity, the central-to-peripheral ratio R^p differs clearly from the ratio (R 1 ^). 

We also compared our NLO and LO calculations of the nuclear modification factor of 
neutral pion production in d+Au collisions, R^,„, in different centrality classes at mid- 
rapidity with the PHENIX data [26] . Within all the given errors in the experimental data, 
the nPDF uncertainties, and the possible differences between the experimental and optical 
Glauber model centrality classes, the EPS09s results are remarkably consistent with the 
centrality systematics. To our knowledge, this is the first time this has been demonstrated. 
Especially, our EPS09s results seem to reproduce the low px slope of the data very well in 
all centrality classes. 

More constraints for the spatial dependence of the nuclear PDFs, and gluons in partic- 
ular, could be obtained from the scheduled p+Pb collisions at the LHC. We demonstrated 
this by calculating the NLO and LO predictions from our framework for the ratio -Rpp b 
in different centrality classes both at mid-rapidity y = and forward rapidity y = 3 for 
y / s/v r /v = 5.0 TeV, which corresponds to the recently achieved p+p cms-energy. 

We believe that the nPDF development presented here is an important step forwards, 
as now a user may for the first time compute the centrality-dependent hard cross-sections 
more consistently with globally analysed nPDFs. Our spatially dependent nPDFs should 
also be applicable in Monte Carlo simulations of nuclear collisions, where the analogues 
of the thickness functions should be straightforwardly obtainable. In addition, our work 
should also give an idea how the future global analyses of nPDFs could be constructed so 
that the spatial dependence would be built in right from the start and not afterwards as 
has been the case here. 

Acknowledgements 

We thank M. Wysocki, J. Rak, T. Lappi, T. Renk and H. Mantysaari for discussions. We 
gratefully acknowledge the following financial support: I.H. from the Magnus Ehrnrooth 
Foundation; K.J.E., I.H. and H.H. from the Academy of Finland, K.J.E.'s Project No. 
133005. C.A.S. is supported by the European Research Council grant HotLHC ERC- 
2001-StG-279579 and by Ministerio de Ciencia e Innovacion of Spain. C.A.S. is a Ramon 
y Cajal researcher. This work (H.H.) was supported in part by the U.S. Department of 
Energy under Grant DE- FG02-93ER40771. 

A Nuclear Collision Geometry 

For clarity, we specify here the modeling and parameters of the nuclear collision geometry, 
i.e. the nuclear thickness functions Ta(s) and Td(s) (see Refs. [57, 58]) and Glauber mod- 
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eling (see Refs. [59, 60]), applied in this study. The calculations of T A (s) and 7a(s) are 
included both in the EKS98s and EPS09s codes. 



A.l Nuclear Thickness Functions 
A. 1.1 Large nuclei 

The total amount of nuclear matter in a colliding nucleus A in the beam direction z at a 
transverse position s is given by the nuclear thickness function 

oo 

T A (s)= J dzp A (s,z), (A.l) 

— oo 

where p(s, z) is the nucleonic number-density of the nucleus, with a normalization conven- 
tion 

A= d 2 sT A (s). (A.2) 

In this study we use the standard two parameter Woods-Saxon density profile for pa, 

n 



PA{s,z) 



1 + exp 



d 



(A.3) 



which is a good approximation for nuclei with A > 4. The parameter values for the 
Woods-Saxon distribution are 

d = 0.54 fm (A.4) 

R A = 1.12A 1 / 3 - 0.86A" 1 / 3 fm, (A.5) 

and for large nuclei the normalization condition (A.2) fixes the constant no as 

"° = ?4(TT(|p- (AJi) 

A. 1.2 Deuterium 

For the thickness function of a deuterium nucleus, the above Woods-Saxon density profile 
is obviously not applicable anymore. Instead, one may formulate this with the deuteron 
wavefunction which describes the probability amplitude for the proton and neutron to 
be separated by a distance r pn . This can be written in terms of the 3 Si- and 3 Di-wave 
components as (see e.g. Ref. [48, 61]) 

<Mr p „) = ^3^(0) + ^3^(0), (A.7) 

where the spin-spherical harmonics yjisityi with S = 1, consist of three components, 

yftiify =(n,m s \LSJM) = Y, (LSM L M s \LSJM)Y L M L {ti)8 ms M s - (A.8) 

Jm s =±l,0 M L ,M S 
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For the radial parts we use the Hulthen form as in [48, 58], 

-/3(ar pn -x c ) 



U 



(r pn ) =NVT 



w(r pn ) =Ne 



1 



1 



-"/(ar pn -x c ) 



e arpn 9(ar pn -x c ) 



(A.9) 
(A.10) 



3(1 - e-T Qr >™) 3(1 _ e -T«-jm)2 
1 + - - + 



or 



pn 



9{ar. 



pn 



where 



2a 



1 — ap 



, in which a 1 = 4.316 fm is related to the experimentally measured 

= 1. For the other 



binding energy, and p is fixed by normalization, f d 3 r pn |^f (r pn )|" 
parameters, obtained by fitting to experimental data, we use the "set 1" quoted in [58]: 

(A.11) 



= 4.680 7 
e = 0.03232 x c 



2.494 




The angular -averaged radial probability distribution for the proton-neutron distance v pn 
in deuteron is given by 



P P n(r P n) = J dn\ifj(r pTl ) 



1 u 2 (r pn ) + w 2 (r pn ) 



-in 



(A.12) 



pn 



For computing the thickness function Td(s) as in Eq. (A.l), we need the nucleon density 
distribution pd( r ) at a distance r from the center of mass of the deuteron. Assuming 
identical proton and neutron masses, we have r = r pn /2. In addition, we require the 
normalization of to be in line with Eq. (A. 2). We thus have 

oo 

T d (s)= J dzp d (s,z), p d (s,z) = WP pn (2r), Jd 2 sT d (s) = 2. (A.13) 

— oo 

A.2 Optical Glauber Model 

Let us then specify the optical Glauber modeling applied for nuclear collisions in this 
study. For further discussion, see e.g. Refs. [59, 60]. Consider first a nucleon-nucleus 
(N+A) collision at an impact parameter b. In the eikonal high collision-energy limit the 
number of binary inelastic collisions is given by 



N» n A (b) = T A (b)a. 



NN 
inel J 



(A.14) 



where Tk(b) is the thickness function defined in Eq. (A.l) and o^Yij is the inelastic nucleon- 
nucleon cross section. One may interpret N^^{h) / A as the probability for an inelastic col- 
lision to take place in the A NN collisions that are possible. Consequently, the probability 
for having no inelastic collisions at all, is 



NN 
el 



(A.15) 
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and the probabililty for at least one inelastic collision becomes 

P^(b) = l-Po(b)-l-e-^ b «. 

The inelastic cross section for the N+A collision we then obtain as 



°£2 = /dV^(b) = |d 2 b(l-e" 



(A.16) 



(A.17) 



As the probability distribution above is expressable in terms of Poissonian probabilities for 
n inelastic collisions, 



= E ^H^>)), (A.18) 

n=l n=l 



1 _ e-^( b ) = V e-^n(b) ( b )] 

71=1 

at a fixed impact parameter we indeed have 

00 



(A.19) 



n=l 



Let us then consider a nucleus-nucleus (j4+-B) collision with collision geometry as 
in Fig. 20. A conveninent choice is to take the impact parameter b along the x axis 
symmetrically around the origin. The transverse density of interacting matter at certain 




A B 

Figure 20. Collision geometry in the transverse plane of the two colliding nuclei. 



impact parameter b can then be computed from the nuclear overlap function, defined as 

T AB (b) = Jd 2 sT A (s 1 )T B (s 2 ), (A.20) 

where si = s + b/2 and S2 = s — b/2. With the normalization for T A (s) in Eq. (A. 2), we 
have 

/ d 2 hT AB {h) = AB. (A.21) 
The number of binary collisions at a given impact parameter b is now 



N^{h)=T AB (h)a^ l 



NN 



(A.22) 
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Analogously to the N+A case above (see e.g. [60]), we may write the probability of an 
inelastic interaction in an A+B collision at an impact parameter b as 

^(b)-l-e^«, (A.23) 

and the inelastic cross section becomes 

^-/^(b) -/d»b(l-.^>. = ). (A,4) 

Figure 21 shows an example of the probability distributions pf n i{h) of Eq. (A.16) for p+Pb 
collisions, and P^fj(b) of Eq. (A.23) for Pb+Pb and d+Pb collisions for cr^?i = 64 mb, 
which corresponds to the cms-energy y / s/v r /v = 2.76 TeV at the LHC. 



90 




|b| [fm] 

Figure 21. The differential inelastic cross section da^Jdb as a function of impact parameter b for 
lead- lead (red solid), dcuteron-lcad (blue long-dashed), and proton- lead (green dashed) collisions 
with = 64 mb. 



The centrality classes in the optical Glauber model can be defined as impact-parameter 
intervals. The "0-c % central" A+B collisions correspond to the most central collisions, 
< b < b c which yield c % of the total inelastic cross section, 

c% = 4b I dV^(b) = °HS; 6c) . (A-25) 
The c\-C2 % centrality class then corresponds to an interval [61,62] for which 

(c 2 -ci)% = 4b[ d 2 b^(b) = U ^ {h AB M) . (A.26) 

01 
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For the studies of different hard-process nuclear modification factors in Sec. 4, we also 
define the average number of binary collisions in a C1-C2 % centrality class: 



(61,62) 



H: tfhT AB {h)a™ 



(7 



AB 
inel 



J? d 2 b 



(A.27) 



where the denominator is simply (02 — c\) % of o - ^. For discussing the analogous centrality 
classes in N+A collisions, we just replace AB by NA in Eqs. (A.25-A.27) above, and also 
T AB by T A in Eq. (A.27). 



B Nuclear modifications r uv and r us 

For completeness, we plot here the nuclear modifications from our fits EPS09sNLOl, 
EPS09sLOl and EKS98s in a lead nucleus for the u valence quarks in Fig. 22 and u 
sea quarks in Fig. 23. The corresponding modifications for gluons are shown in Fig. 5. 





Figure 22. The spatially dependent modification of the u distribution in a lead nucleus, 
rP b (.T, Q 2 ,s), from EPS09sNLOl (upper left), EPS09sLOl (upper right) and EKS98s(lower plot) 
as a function of x and s at the initial scale Q 2 = 1.69(2.25) GcV 2 of EPS09 (EKS98). 
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Figure 23. The spatially dependent modification of the uy distribution in a lead nucleus, 
r F g h (x,Q 2 ,s), from EPS09sNLOl (upper left), EPS09sLOl (upper right) and EKS98s(lower plot) 
as a function of x and s at the initial scale Q 2 = 1.69 (2.25) GcV 2 of EPS09 (EKS98). 
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